home *** CD-ROM | disk | FTP | other *** search
/ Software of the Month Club 1997 February / Software of the Month Club 1997 February.iso / pc / dos / edu / crystal / crystal6.cpp next >
Encoding:
C/C++ Source or Header  |  1996-06-28  |  12.0 KB  |  330 lines
  1. // crystal.cpp
  2. // A program to calculate spacing, angles between planes (hkl)
  3. // and angles between zones [uvw]
  4. // copyright May 1996, by Gordon Vrdoljak
  5. #include <iostream.h>
  6. #include <math.h>
  7.  
  8. double spacing(double a, double b, double c, double alpha, double beta,
  9.         double gamma, int h, int k, int l, char system);
  10.  
  11. double planar_angle(double a, double b, double c, double alpha, double beta,
  12.         double gamma, int h1, int k1, int l1, int h2, int k2, int l2,
  13.         char system );
  14. double zonal_angle(double a, double b, double c, double alpha, double beta,
  15.         double gamma, int u1, int v1, int w1, int u2, int v2, int w2,
  16.         char system );
  17.  
  18. main()
  19. {
  20. // variable declarations:
  21. char system, query;
  22. double a, b, c, alpha, beta, gamma, d, phi, rho;
  23. int choice, h, k, l, h1, k1, l1;
  24. int h2, k2, l2, u1, v1, w1, u2, v2, w2;
  25.  
  26. cout << "A program to calculate the spacing between the" << endl
  27.         << "crystallographic (hkl) planes, angle between planes," << endl
  28.         << "and the angle between zones.\n" << endl;
  29. cout << "Is the material:\n" << "\n\t(1) Cubic\n"
  30.         << "\t(2) Tetragonal\n" << "\t(3) Orthorhomic\n"
  31.         << "\t(4) Hexagonal\n" << "\t(5) Rhombohedral\n"
  32.         << "\t(6) Monoclinic\n" << "\t(7) Triclinic";
  33. cout << "\nInput number (1-7):";
  34. cin >> system;
  35. switch(system)
  36.     {
  37.     case '1':
  38.         cout << "\nInput value for 'a' cell parameter:";
  39.         cin >> a;
  40.         b = a; c = a; alpha = 1.570796; beta = alpha; gamma = alpha;
  41.         break;
  42.     case '2':
  43.         cout << "\nInput value for 'a' cell parameter:";
  44.         cin >> a;
  45.         cout << "\nInput value for 'c' cell parameter:";
  46.         cin >> c;
  47.         b = a; alpha = 1.570796; beta = alpha; gamma = alpha;
  48.         break;
  49.     case '3':
  50.         cout << "\nInput value for 'a' cell parameter:";
  51.         cin >> a;
  52.         cout << "\nInput value for 'b' cell parameter:";
  53.         cin >> b;
  54.         cout << "\nInput value for 'c' cell parameter:";
  55.         cin >> c;
  56.         alpha = 1.570796; beta = alpha; gamma = alpha;
  57.         break;
  58.     case '4':
  59.         cout << "\nInput value for 'a' cell parameter:";
  60.         cin >> a;
  61.         cout << "\nInput value for 'c' cell parameter:";
  62.         cin >> c;
  63.         b = a; alpha = 1.570796; beta = alpha; gamma = 2.094395;
  64.         break;
  65.     case '5':
  66.         cout << "\nInput value for 'a' cell parameter:";
  67.         cin >> a;
  68.         cout << "\nInput value for 'alpha' in degrees:";
  69.         cin >> alpha;
  70.         alpha = alpha * 3.141593 / 180;
  71.         b = a; c = a; beta = alpha; gamma = alpha;
  72.         break;
  73.     case '6':
  74.         cout << "\nInput value for 'a' cell parameter:";
  75.         cin >> a;
  76.         cout << "\nInput value for 'b' cell parameter:";
  77.         cin >> b;
  78.         cout << "\nInput value for 'c' cell parameter:";
  79.         cin >> c;
  80.         cout << "\nInput value for 'beta' in degrees:";
  81.         cin >> beta;
  82.         beta = beta * 3.141593 / 180;
  83.         alpha = 1.570796; gamma = 1.570796;
  84.         break;
  85.     case '7':
  86.         cout << "\nInput value for 'a' cell parameter:";
  87.         cin >> a;
  88.         cout << "\nInput value for 'b' cell parameter:";
  89.         cin >> b;
  90.         cout << "\nInput value for 'c' cell parameter:";
  91.         cin >> c;
  92.         cout << "\nInput value for 'alpha' in degrees:";
  93.         cin >> alpha;
  94.         alpha = alpha * 3.141593 / 180;
  95.         cout << "\nInput value for 'beta' in degrees:";
  96.         cin >> beta;
  97.         beta = beta * 3.141593 / 180;
  98.         cout << "\nInput value for 'gamma' in degrees:";
  99.         cin >> gamma;
  100.         gamma = gamma * 3.141593 / 180;
  101.         break;
  102.         }
  103. cout << "\nDo you want to find:" << "\n\t (1) the inter planar spacing, d"
  104.         << "\n\t (2) the angle between planes (h1 k1 l1) \n\t\t and (h2 l2 k2), or"
  105.         << "\n\t (3) the interzonal angle between zones \n\t\t (u1 v1 w1) and"
  106.         << "(u2 v2 w2)"
  107.         << "\n\nInput (1, 2, 3): ";
  108. cin >> choice;
  109. switch(choice)
  110.     {
  111.     case 1:
  112.         cout << "\nInput h, k, l of plane." << "\nh:";
  113.         cin >> h; cout << "\nk:"; cin >> k; cout << "\nl:"; cin >> l;
  114.         d = spacing(a, b, c, alpha, beta, gamma, h, k, l, system);
  115.         cout << "\nthe d spacing is: " << d << " Angstroms\n";
  116.         break;
  117.     case 2:
  118.         if (system == '5')
  119.             cout << " Please convert to hexagonal indices and use "
  120.                   << "the hexagonal system\n for this calculation.\n";
  121.         else {
  122.             cout << "\nInput h, k, l of plane 1." << "\nh1:"; cin >> h1;
  123.             cout << "\nk1:"; cin >> k1; cout << "\nl1:"; cin >> l1;
  124.             cout << "\nInput h, k, l of plane 2." << "\nh2:"; cin >> h2;
  125.             cout << "\nk2:"; cin >> k2; cout << "\nl2:"; cin >> l2;
  126.             phi = planar_angle(a, b, c, alpha, beta, gamma, h1, k1, l1, h2, k2, l2,
  127.                 system);
  128.             cout << "\nThe angle between " << h1 << k1 << l1 << " and " << h2 << k2
  129.                 << l2 << " is: " << phi << "radians,\n or " << 180 * phi /
  130.                 3.141592654    << " degrees.\n";
  131.             }
  132.         break;
  133.     case 3:
  134.         if (system == '5')
  135.             cout << " Please convert to hexagonal indices and use "
  136.                   << "the hexagonal system\n for this calculation.\n";
  137.         else {
  138.             cout << "\nInput u, v, w of zone 1." << "\nu1:"; cin >> u1;
  139.             cout << "\nv1:"; cin >> v1; cout << "\nw1:"; cin >> w1;
  140.             cout << "\nInput u, v, w of plane 2." << "\nu2:"; cin >> u2;
  141.             cout << "\nv2:"; cin >> v2; cout << "\nw2:"; cin >> w2;
  142.             rho = zonal_angle(a, b, c, alpha, beta, gamma, u1, v1, w1, u2, v2, w2,
  143.                 system);
  144.             cout << "\nThe angle between " << u1 << v1 << w1 << " and " << u2 << v2
  145.                 << w2 << " is: " << rho << "radians,\n or " << 180 * rho / 3.141592654
  146.                 << " degrees.\n";
  147.         }
  148.         break;
  149.     }
  150. cout << "a = " << a << ", b = " << b << ", c = " << c;
  151. cout << "\nalpha = " << alpha << ", beta = " << beta
  152.         << ", gamma = " << gamma;
  153. return 0;
  154. }
  155.  
  156. // Function to calculate d spacing
  157. double spacing(double a, double b, double c, double alpha, double beta,
  158.         double gamma, int h, int k, int l, char system)
  159. {
  160.     double d, q;
  161.     switch(system)
  162.         {
  163.         case '1':
  164.             d = sqrt( pow (a,  2) / (h * h + k * k + l * l));
  165.             return d;
  166.         case '2':
  167.             q = (h * h + k * k) / (a * a) + l * l / (c * c);
  168.             d = sqrt( 1 / q);
  169.             return d;
  170.         case '3':
  171.             q = h * h / a * a + k * k / b * b + l * l / c * c;
  172.             d = sqrt( 1 / q);
  173.             return d;
  174.         case '4':
  175.             q = 4 * (h * h + h * k + k * k) / 3 * a * a + l * l / c * c;
  176.             d = sqrt( 1 / q);
  177.             return d;
  178.         case '5':
  179.             q = (1 / a * a)*((1 + cos(alpha)) * ((h * h + k * k + l * l) -
  180.                 (1 - pow(tan (0.5 * alpha),2) ) * (h * k + k * l + l * h)) / (1 +
  181.                 cos (alpha) - 2 * pow(cos (alpha), 2)));
  182.             d = sqrt( 1 / q);
  183.             return d;
  184.         case '6':
  185.             q = (1 / a * a) * h * h / (pow(sin(beta),2)) + k * k / (b * b) +
  186.                 l * l / (c * c * pow(sin(beta), 2)) - 2 * h * l * cos (beta) /
  187.                 (a * c * pow(sin(beta), 2));
  188.             d = sqrt( 1 / q);
  189.             return d;
  190.         case '7':
  191.             q = (1 / (a * a * b * b * c * c * (1 - pow(cos(alpha), 2) -
  192.                 pow(cos(beta), 2) - pow(cos(gamma), 2) + 2 * cos(alpha) * cos(beta)
  193.                  * cos(gamma)))) * (b * b * c * c * pow(sin(alpha), 2) * h * h
  194.                  + a * a * c * c * pow(sin(beta), 2) * k * k + a * a * b * b *
  195.                  pow(sin(gamma), 2) * l * l + 2 * a * b * c * c * (cos(alpha) *
  196.                  cos(beta) - cos(gamma)) * h * k + 2 * a * a * b * c * (cos(beta)
  197.                  * cos(gamma) - cos(alpha)) * k * l + 2 * a * b * b * c *
  198.                  (cos(gamma) * cos(alpha) - cos(beta)) * l * h);
  199.             d = sqrt( 1 / q);
  200.             return d;
  201.         }
  202. }
  203.  
  204. // Function to calculate interplanar angles
  205. double planar_angle(double a, double b, double c, double alpha, double beta,
  206.     double gamma, int h1, int k1, int l1, int h2, int k2, int l2, char system )
  207. {
  208. double phi;
  209.     switch(system)
  210.         {
  211.         case '1':
  212.             phi = acos((h1 * h2 + k1 * k2 + l1 * l2) / sqrt((h1 * h1 + k1 * k1 +
  213.                         l1 * l1) * (h2 * h2 + k2 * k2 + l2 * l2)));
  214.             return phi;
  215.         case '2':
  216.  
  217.             phi = acos(((h1 * h2 + k1 * k2) / (a * a) + l1 * l2 / (c * c)) /
  218.                         sqrt((h1 * h1 + k1 * k1) / (a * a) + l1 * l1 / (c * c)) * ((
  219.                         h2 * h2 + k2 * k2) / (a * a) + l2 * l2 / (c * c)));
  220.             return phi;
  221.         case '3':
  222. //            double phi;
  223.             phi = acos((h1 * h2 / (a * a) + k1 * k2 / (b * b) + l1 * l2 / (c * c))
  224.                         / sqrt((h1 * h1 / (a * a) + k1 * k1 / (b * b) + l1 * l1 /
  225.                         (c * c)) * (h2 * h2 / (a * a) + k2 * k2 / (b * b) + l2 * l2 /
  226.                         (c * c))));
  227.             return phi;
  228.         case '4':
  229. //            double phi;
  230.             phi = acos((h1 * h2 + k1 * k2 + 0.5 * (h1 * k2 + k1 * h2) + 0.75 * a
  231.                         * a * l1 * l2 / (c * c)) / sqrt((h1 * h1 + k1 * k1 + h1 * k1
  232.                         + 0.75 * a * a * l1 * l1 / (c * c)) * (h2 * h2 + k2 * k2 + h2
  233.                         * k2 + 0.75 * a * a * l2 * l2 / (c * c))));
  234.             return phi;
  235.         case '5':
  236. //            double phi;
  237.             phi = 999;
  238.             return phi;
  239.         case '6':
  240. //            double phi;
  241.             phi = acos((h1 * h2 / (a * a) + k1 * k2 / (b * b) * pow(sin(beta), 2)
  242.                         + l1 * l2 / (c * c) - (l1 * h2 + l2 * h1) * cos(beta)) / sqrt(
  243.                         (h1 * h1 / (a * a) + k1 * k1 / (b * b) * pow(sin(beta), 2) +
  244.                         l1 * l1 / (c * c) - 2 * h1 *l1 * cos(beta) / (a * c)) * (h2 *
  245.                         h2 / (a * a) + k2 * k2 / (b * b) * pow(sin(beta), 2) + l2 *
  246.                         l2 / (c * c) - 2 * h2 * l2 / (a * c) * cos(beta))));
  247.             return phi;
  248.         case '7':
  249. //            double phi;
  250.             phi = acos((h1 * h2 * b * b * c * c * pow(sin(alpha), 2) + k1 * k2 * a
  251.                         * a * c * c * pow(sin(beta), 2) + l1 * l2 * a * a * b * b *
  252.                         pow(sin(gamma), 2) +  a * b * c * c * (cos(alpha) * cos(beta)
  253.                         - cos(gamma)) * (k1 * h2 + h1 * k2) + a * b * b * c * (
  254.                         cos(gamma) * cos(alpha) - cos(beta)) * (h1 * l2 + l1 * h2) +
  255.                         a * a * b * c * (cos(beta) * cos(gamma) - cos(alpha)) * (k1 *
  256.                         l2 + l1 * k2)) / (
  257.                         sqrt(h1 * h1 * b * b * c * c * pow(sin(alpha), 2) + k1 * k1 *
  258.                         a * a * c * c * pow(sin(beta), 2) + l1 * l1 * a * a * b * b *
  259.                         pow(sin(gamma), 2) + 2 * h1 * k1 * a * b * c * c * (cos(alpha)
  260.                         * cos(beta) - cos(gamma)) + 2 * h1 * l1 * a * b * b * c * (
  261.                         cos(gamma) * cos(alpha) - cos(beta)) + 2 * k1 * l1 * a * a *
  262.                         b * c * (cos(beta) * cos(gamma) - cos(alpha))) *
  263.                         sqrt(h2 * h2 * b * b * c * c * pow(sin(alpha), 2) + k2 * k2 *
  264.                         a * a * c * c * pow(sin(beta), 2) + l2 * l2 * a * a * b * b *
  265.                         pow(sin(gamma), 2) + 2 * h2 * k2 * a * b * c * c * (cos(alpha)
  266.                         * cos(beta) - cos(gamma)) + 2 * h2 * l2 * a * b * b * c * (
  267.                         cos(gamma) * cos(alpha) - cos(beta)) + 2 * k2 * l2 * a * a *
  268.                         b * c * (cos(beta) * cos(gamma) - cos(alpha)))
  269.                         ));
  270.             return phi;
  271.         }
  272. }
  273.  
  274.  
  275. //      float zonal_angle(float a, float b, float c, float alpha, float beta, float gamma,
  276. //        int u1, int v1, int w1, int u2, int v2, int w2, char system );
  277. // Function to calculate interzonal angles
  278. double zonal_angle(double a, double b, double c, double alpha, double beta,
  279.         double gamma, int u1, int v1, int w1, int u2, int v2, int w2,
  280.         char system )
  281. {
  282. double rho;
  283.     switch(system)
  284.         {
  285.         case '1':
  286.             rho = acos((u1 * u2 + v1 * v2 + w1 *w2) /
  287.                     sqrt((u1 * u1 + v1 * v1 + w1 * w1) * (u2 * u2 + v2 * v2 + w2 *
  288.                     w2)));
  289.             return rho;
  290.         case '2':
  291.             rho = acos((a * a * (u1 * u2 + v1 * v2) + c * c * w1 * w2) /
  292.                     sqrt((a * a * (u1 * u1 + v1 * v1) + c * c * w1 * w1) *
  293.                     (a * a * (u2 * u2 + v2 * v2) + c * c * w2 * w2)));
  294.             return rho;
  295.         case '3':
  296.             rho = acos((a * a * u1 * u2 + b * b * v1 * v2 + c * c * w1 * w2) /
  297.                     sqrt((a * a * u1 * u1 + b * b * v1 * v1 + c * c * w1 * w1) *
  298.                     (a * a * u2 * u2 + b * b * v2 * v2 + c * c * w2 * w2)));
  299.             return rho;
  300.         case '4':
  301.             rho = acos((u1 * u2 + v1 * v2 - 0.5 * (u1 * v2 + v1 * u2) + c * c / (
  302.                     a * a) * w1 * w2) /
  303.                     sqrt((u1 * u1 + v1 * v1 - u1 * v1 + c * c / ( a * a) * w1 * w1) *
  304.                     (u2 * u2 + v2 * v2 - u2 * v2 + c * c / (a * a) * w2 * w2)));
  305.             return rho;
  306.         case '5':
  307.             rho = 999;
  308.             return rho;
  309.         case '6':
  310.             rho = acos((a * a * u1 * u2 + b * b * v1 * v2 + c * c * w1 * w2 + a *
  311.                     c * (w1 * u2 + u1 * w2) * cos(beta)) /
  312.                     sqrt((a * a * u1 * u1 + b * b * v1 * v1 + c * c * w1 * w1 + 2 * a
  313.                     * c * u1 * w1 * cos(beta)) *
  314.                     (a * a * u2 * u2 + b * b * v2 * v2 + c * c * w2 * w2 + 2 * a * c
  315.                     * u2 * w2 * cos(beta))));
  316.             return rho;
  317.         case '7':
  318.             rho = acos((a * a * u1 * u2 + b * b * v1 * v2 + c * c * w1 * w2 + b *
  319.                     c * (v1 * w2 + w1 * v2) * cos(alpha) + a * c * (w1 * u2 + u1 *
  320.                     w2) * cos(beta) + a * b * (u1 * v2 + v1 * u2) * cos(gamma)) /
  321.                     (sqrt(a * a * u1 * u1 + b * b * v1 * v1 + c * c * w1 * w1 + 2 *
  322.                     b * c * v1 * w1 * cos(alpha) + 2 * c * a * w1 * u1 * cos(beta) +
  323.                     2 * a * b * u1 * v1 * cos(gamma)) *
  324.                     sqrt(a * a * u2 * u2 + b * b * v2 * v2 + c * c * w2 * w2 + 2 *
  325.                     b * c * v2 * w2 * cos(alpha) + 2 * c * a * w2 * u2 * cos(beta) +
  326.                     2 * a * b * u2 * v2 * cos(gamma))));
  327.         return rho;
  328.         }
  329. }
  330.